sc10x analysis

load dependancies

Baf53cre Ahr-cko vs ctl,
CR infection,
Intestinal RORgt+ Immune Cells

loading 72,000 for all
cellranger calling 34,620 cells, mean reads 11,163

t1 advanced

load obj

GEX.seur <- readRDS("./20230927_10x_SZJ.new20270809.preAnno.rds")
GEX.seur
## An object of class Seurat 
## 16621 features across 24449 samples within 1 assay 
## Active assay: RNA (16621 features, 1224 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
#
color.FB <- ggsci::pal_d3("category20c")(20)[c(2,7,12,17,
                                               1,6,11,16
                                               )]
  
color.cnt <- color.FB[c(5,1)]
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno", label.size = 2.75) +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

filtering

# remove Mix
# remove DF0.1-Doublets except cycling ILC3/Th
GEX.seur <- subset(GEX.seur, subset = preAnno != "Mix")
GEX.seur <- subset(GEX.seur, subset = DoubletFinder0.1=="Singlet" | preAnno %in% c("Th.CC","ILC3.CC"))
GEX.seur
## An object of class Seurat 
## 16621 features across 22088 samples within 1 assay 
## Active assay: RNA (16621 features, 1224 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno", label.size = 2.75) +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

VlnPlot(subset(GEX.seur, subset = cnt %in% c("CTL")), 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "preAnno")

VlnPlot(subset(GEX.seur, subset = cnt %in% c("CKO")), 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "preAnno")

markers.new <-     c("Cd3d","Cd3e","Cd4","Cd8a",
                     "Klf2","Ccr7","Lef1","S1pr1",
                     "Tubb2a","Rflnb","Igfbp4","Sell",
                     "Itm2a","Itga6","Smc4",#"Sox4","Pdcd1",
                     "Ptprc",#"Fam241a","Slamf6",
                     "Izumo1r","Tbc1d4",
                     "Ctla4","Il10", 
                     "Hopx","Foxp3","Ccr2","Maf",
                     "Capg","Slamf1","Icos",
                     "Il17a","Nebl","S100a11","S100a10",
                     "Top2a","Mki67","Pcna","Mcm6",
                     
                     "Fasl","Gzma","Gzmk","Ccr5",
                     "Il12rb2",
                     "Jun","Fos","Ier2","Egr1",
                     
                     "Ifng", "Tnf",
                     "Tbx21",
                     "Klrb1c","Xcl1","Nkg7","Plac8","Cd160",
                     
                     
                     
                     "Ccl5",
                     "Klrd1","Cd7","Itgae",
                     
                     "Tcrg-C1","Trdv4","Cd163l1","Blk","Mmp25",
                     
                     "Calca","Gata3","Il17rb","Il13",
                     "Csf2","Dgat2",
                     "Il4","Hilpda","Ar","Il1rl1",
                     "Il5",#"Cd24a",
                     "Ccl1","Areg","Ahr",
                     "Cxcl10",
                     "Gbp6","Stat1","Isg15","Ifit1",
                     
                     "Ncr1","Gzmc","Gzmb","Irf8",
                     "Cxcr6","Irf7","Slc16a6","Klrk1",
                     "Tmem176b","Fcer1g","Car2","Ckb",
                     "Cd69","Jaml","Cxcl9",
                     "Il22","Gem","Pxdc1","Cdc14a",
                     "Sdc4","Vegfa","Bst2","Nmrk1",
                     "Hmgn3","Cd81","Cx3cl1",
                     "Cd74","H2-Aa","H2-Eb1","H2-Ab1",
                     
                     
                     "Ccr6","Batf3","Rorc","Maff",
                     "Icam1","Il18r1",
                     
                     "Il17f"
                     
                     )

DotPlot(GEX.seur, features = markers.new, group.by = "preAnno",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 9)) + 
  scale_y_discrete(limits=rev)

re-clustering

GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1500)

# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2

head(VariableFeatures(GEX.seur), 200)
##   [1] "Hist1h1b"  "Il22"      "Ccl4"      "Il17a"     "Gzma"      "Ccl1"     
##   [7] "Ccl5"      "Cxcl2"     "Cxcl10"    "Areg"      "Xcl1"      "Hist1h2ae"
##  [13] "Calca"     "Stmn1"     "Cxcl1"     "Hspa1a"    "S100g"     "Ube2c"    
##  [19] "Il13"      "Il5"       "Ccl20"     "Gzmc"      "Il4"       "Cd74"     
##  [25] "Pclaf"     "S100a6"    "Gzmb"      "Cxcl3"     "Cxcl9"     "Tnfrsf4"  
##  [31] "Il17f"     "Ifng"      "Lgals1"    "Ccl3"      "Ifi27l2a"  "Hs3st1"   
##  [37] "Il10"      "Trdc"      "Egr1"      "Penk"      "H2-Aa"     "Top2a"    
##  [43] "Mki67"     "Stmn2"     "Csf2"      "Il2"       "Pcp4"      "Hist1h2ap"
##  [49] "Tph1"      "Birc5"     "Plac8"     "Defa24"    "H2-Eb1"    "Hmgb2"    
##  [55] "Hist1h3c"  "Cenpf"     "Serpine1"  "Sostdc1"   "Hspb1"     "Ctla2a"   
##  [61] "Jun"       "Egr2"      "Hist1h2ab" "Serpinb1a" "Lgals3"    "Ifitm2"   
##  [67] "Rrm2"      "Try5"      "S100a4"    "Ifitm1"    "Bambi"     "Cd36"     
##  [73] "Tmsb4x"    "Rrad"      "Ccr7"      "Ctla4"     "Actb"      "Crip1"    
##  [79] "Tubb5"     "Hspa1b"    "Ifit1"     "Ly6a"      "Hist1h1e"  "Gzmf"     
##  [85] "Bcl2a1b"   "Cdca8"     "Fabp5"     "Gzme"      "Nusap1"    "H2-Ab1"   
##  [91] "Tubb2a"    "Akap12"    "Klf2"      "Vim"       "G0s2"      "Mt1"      
##  [97] "Dppa2"     "Nkg7"      "Rgcc"      "Defa21"    "Tnfsf8"    "Igkc"     
## [103] "Fos"       "Apoe"      "Spp1"      "C2cd4b"    "Cd7"       "Ptgs2"    
## [109] "Tyrobp"    "Egr3"      "Ifitm3"    "Lyz1"      "Hbegf"     "Igha"     
## [115] "Ccnb2"     "Dnaja1"    "Cldn5"     "Hes1"      "Cd83"      "Gadd45g"  
## [121] "Ms4a4b"    "Cks1b"     "Gata3"     "Actg2"     "Rasl11a"   "Gm10827"  
## [127] "Defa22"    "Gadd45b"   "Isg15"     "Gfod2"     "Tyms"      "Clspn"    
## [133] "Jchain"    "Trac"      "Gzmd"      "Trbc1"     "Klra5"     "Iigp1"    
## [139] "Hopx"      "Tuba1b"    "Ptpn13"    "Arg1"      "Hmmr"      "Hba-a1"   
## [145] "Il21"      "Defa30"    "Ermn"      "Fbxo5"     "Ccna2"     "Dennd5b"  
## [151] "Hilpda"    "Klra1"     "Muc3"      "Atf3"      "Spc24"     "Odc1"     
## [157] "Klrd1"     "AY761184"  "Il1r2"     "Slpi"      "Dnajb1"    "Il1rl1"   
## [163] "Cdca3"     "Acod1"     "Cd3g"      "Fst"       "Tpx2"      "Klra7"    
## [169] "H2afz"     "Dut"       "Bcl2a1a"   "Kif11"     "Hist1h1d"  "Cd5"      
## [175] "Il6"       "Cd24a"     "Zfp36"     "Gm4841"    "Serpinb9"  "Igfbp4"   
## [181] "Fabp4"     "Cd163l1"   "Trdv4"     "Esco2"     "Socs2"     "Cenpe"    
## [187] "Egr4"      "Chad"      "Eprs"      "Gm14851"   "Serpinb6b" "Smc2"     
## [193] "Id2"       "Izumo1r"   "Arl5c"     "Gem"       "Tent5a"    "Ccr2"     
## [199] "Akr1c18"   "Cd3d"
GEX.seur <- ScaleData(GEX.seur, features = rownames(GEX.seur))
## Centering and scaling data matrix
# exclude MT genes  and more 

DIG <- grep("^mt-|^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^AC|^AA|^AW|^BC|^Gm|^Hist|Rik$|-ps|^Ifi|^Isg|^Egr",    
            # try adding some few to reduce the mixed Ifit+ 
            rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))


VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
                                                c(#MT_gene,
                                                  DIG, 
                                                  CC_gene) )

GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur))
## PC_ 1 
## Positive:  Fcer1g, Tmem176a, Gem, Tmem176b, Itm2b, Pxdc1, Sdc4, Il22, Car2, Nrgn 
##     Gpr183, Cdkn1a, Eprs, Zfp36l1, Ckb, Prr29, Ccdc184, Klrb1b, Hk2, Odc1 
##     Klf4, Cwc25, Tnfsf11, Rgs1, Atf3, Fam110a, Arrdc4, Gpx1, Dad1, Ostf1 
## Negative:  Cd3d, Cd3g, Cd2, Gramd3, Ms4a4b, Fyb, Gimap6, Cd28, Cd3e, Ms4a6b 
##     Il21r, Dusp10, Lef1, Klf2, Lat, Ccr7, Pmepa1, Slfn2, Satb1, S1pr1 
##     Stk4, Gimap4, Gimap1, Cd5, Cd27, Coro1a, Cd247, Sh2d2a, Got1, Gimap7 
## PC_ 2 
## Positive:  Ccr7, Vegfa, Cdc14a, Hmgn3, Cd81, Bst2, S1pr1, Sdc4, Ramp1, Cx3cl1 
##     Klf2, Cd74, Batf3, Igfbp4, Lef1, Cryaa, Rflnb, Bmp2, Irs2, Rgcc 
##     Myof, Nrp1, Ccr6, Cited2, Ttn, H2-Aa, Lmna, Cd83, Il22, Tgm2 
## Negative:  Tmsb4x, Bcl2a1d, Cxcr6, Serpinb9, Pfn1, Lgals1, Capg, Serpinb6b, Ucp2, Pik3r1 
##     Icos, Ctsw, Cd52, Serpinb1a, Bcl2a1b, Actn2, Srgn, Ltb4r1, Slc16a6, Ikzf2 
##     Ptprcap, Maf, Clic1, Slc7a6os, S100a6, Rgs1, Nkg7, S100a11, Rbpms2, Arl5c 
## PC_ 3 
## Positive:  Hs3st1, Calca, Ptpn13, Il5, Gata3, Ccl1, Klrg1, Il17rb, Il13, Il1rl1 
##     Ipmk, Csf2, Areg, Il4, Slc7a8, Otulinl, Tspan13, Deptor, Lmo4, Klf5 
##     Tnfrsf18, Hes1, Spcs3, Rab27b, Tent5a, Mras, Rab4a, Ptgir, Ar, Dgat2 
## Negative:  Fcer1g, Ckb, Ms4a4b, Lef1, Car2, Ccr7, Serpinb1a, Gpx1, Ms4a6b, Gimap4 
##     S1pr1, Klf2, Tubb2a, Klrk1, Gramd3, Capg, Cd3d, Serpinb6b, Serpinb9, Rflnb 
##     Prr29, Stmn1, Klrb1b, Il22, Upp1, Calm1, Rbpms2, Cd2, Rgcc, Slc7a6os 
## PC_ 4 
## Positive:  Pclaf, Spc24, Actb, Stmn1, Tnfrsf4, Ccna2, Asf1b, Tk1, Racgap1, Kif4 
##     Vim, S100a4, S100a6, Ccr6, Ctla4, Knl1, Smc2, Actg1, Lgals1, Bst2 
##     Mcm3, Cenpm, Lig1, Hmgn3, Cd74, Dut, Ramp1, Il17a, Ly6a, Tubb5 
## Negative:  Serpinb9, Serpinb6b, Pik3r1, Ctsw, Lef1, Upp1, Klrk1, Satb1, Slc7a6os, Rbpms2 
##     Tubb2a, Ncr1, Rflnb, Gimap6, Gimap4, Irf7, Dusp10, Ikzf2, Pcp4, Slc16a10 
##     Slc16a6, Serpinb1a, Klf2, Sell, Lmo4, Ccr7, Gjb2, Ppfia1, Cd200r4, Trf 
## PC_ 5 
## Positive:  Bcl2a1b, Cd3g, Ccr2, Ctla4, Cd3e, Hopx, Ccl5, S100a6, Tnfrsf4, S100a10 
##     Ccr5, Got1, Fth1, Cd6, Actb, Il10, Odc1, Cd52, Actn2, Cd163l1 
##     S100a11, Cd3d, Foxp3, Sh2d2a, Plac8, Ltb, Blk, Ubald2, Nkg7, Cd40lg 
## Negative:  Pclaf, Stmn1, Spc24, Ccna2, Asf1b, Kif4, Knl1, Cenpm, Tk1, Smc2 
##     Shcbp1, Spc25, Tubb5, Lef1, Racgap1, Ccr7, Nrm, Kif15, S1pr1, Rflnb 
##     Cdkn3, Tubb2a, Pbk, Tuba1b, Cenph, Aspm, Mcm3, Hmgn2, Esco2, Hs3st1
length(setdiff(VariableFeatures(object = GEX.seur),
                                                c(#MT_gene,
                                                  DIG,CC_gene)))
## [1] 1197
setdiff(VariableFeatures(object = GEX.seur),
                                                c(#MT_gene,
                                                  DIG,CC_gene))[1:300]
##   [1] "Il22"      "Ccl4"      "Il17a"     "Gzma"      "Ccl1"      "Ccl5"     
##   [7] "Cxcl2"     "Cxcl10"    "Areg"      "Xcl1"      "Calca"     "Stmn1"    
##  [13] "Cxcl1"     "S100g"     "Il13"      "Il5"       "Ccl20"     "Gzmc"     
##  [19] "Il4"       "Cd74"      "Pclaf"     "S100a6"    "Gzmb"      "Cxcl3"    
##  [25] "Cxcl9"     "Tnfrsf4"   "Il17f"     "Ifng"      "Lgals1"    "Ccl3"     
##  [31] "Hs3st1"    "Il10"      "Penk"      "H2-Aa"     "Stmn2"     "Csf2"     
##  [37] "Il2"       "Pcp4"      "Tph1"      "Plac8"     "Defa24"    "H2-Eb1"   
##  [43] "Serpine1"  "Sostdc1"   "Ctla2a"    "Serpinb1a" "Lgals3"    "Try5"     
##  [49] "S100a4"    "Bambi"     "Cd36"      "Tmsb4x"    "Rrad"      "Ccr7"     
##  [55] "Ctla4"     "Actb"      "Crip1"     "Tubb5"     "Ly6a"      "Gzmf"     
##  [61] "Bcl2a1b"   "Fabp5"     "Gzme"      "H2-Ab1"    "Tubb2a"    "Akap12"   
##  [67] "Klf2"      "Vim"       "G0s2"      "Mt1"       "Dppa2"     "Nkg7"     
##  [73] "Rgcc"      "Defa21"    "Tnfsf8"    "Apoe"      "Spp1"      "C2cd4b"   
##  [79] "Cd7"       "Ptgs2"     "Tyrobp"    "Lyz1"      "Hbegf"     "Cldn5"    
##  [85] "Hes1"      "Cd83"      "Gadd45g"   "Ms4a4b"    "Gata3"     "Actg2"    
##  [91] "Rasl11a"   "Defa22"    "Gadd45b"   "Gfod2"     "Gzmd"      "Klra5"    
##  [97] "Iigp1"     "Hopx"      "Tuba1b"    "Ptpn13"    "Arg1"      "Il21"     
## [103] "Defa30"    "Ermn"      "Fbxo5"     "Ccna2"     "Dennd5b"   "Hilpda"   
## [109] "Klra1"     "Muc3"      "Atf3"      "Spc24"     "Odc1"      "Klrd1"    
## [115] "Il1r2"     "Slpi"      "Il1rl1"    "Acod1"     "Cd3g"      "Fst"      
## [121] "Klra7"     "H2afz"     "Dut"       "Bcl2a1a"   "Cd5"       "Il6"      
## [127] "Cd24a"     "Zfp36"     "Serpinb9"  "Igfbp4"    "Fabp4"     "Cd163l1"  
## [133] "Esco2"     "Socs2"     "Chad"      "Eprs"      "Serpinb6b" "Smc2"     
## [139] "Id2"       "Izumo1r"   "Arl5c"     "Gem"       "Tent5a"    "Ccr2"     
## [145] "Akr1c18"   "Cd3d"      "Lef1"      "Sptssb"    "Gramd3"    "Blnk"     
## [151] "Nrgn"      "Lyz2"      "Tnfsf11"   "Lgals7"    "Klra6"     "Lmo4"     
## [157] "Ier3"      "Nlrp3"     "Adra2a"    "Bcl2a1d"   "Cd2"       "Upp1"     
## [163] "Vegfa"     "Cd52"      "Adm"       "Prss2"     "Tpm2"      "Cd28"     
## [169] "Tnfrsf9"   "Dusp14"    "Basp1"     "Rsad2"     "Serpine2"  "Rflnb"    
## [175] "Actn2"     "Gcg"       "Msx1"      "Mcm3"      "Ptprz1"    "Cst8"     
## [181] "Il1rn"     "Evx1os"    "Nkx6-2"    "Rgs1"      "Pcdh10"    "Zg16"     
## [187] "Olfm4"     "Ptcra"     "Atp1b2"    "Hmgn2"     "Tulp3"     "Gpr83"    
## [193] "Klre1"     "Cdkn1a"    "Irf7"      "Cebpd"     "Timp3"     "S1pr1"    
## [199] "Lmna"      "Peg3"      "Ccr3"      "Emp3"      "Rnd3"      "Aldh1a3"  
## [205] "Foxf1"     "H1f0"      "Rgs16"     "Gla"       "Pfn1"      "Kif4"     
## [211] "Sox4"      "Nmrk1"     "Palld"     "Coro1a"    "Aspm"      "H2afx"    
## [217] "Kcnq1ot1"  "Dusp2"     "Plk2"      "Plek"      "Asb4"      "Bst2"     
## [223] "Gbp2"      "Hoxa7"     "Tshz2"     "Ccr5"      "Asf1b"     "Clu"      
## [229] "Foxd1"     "Fzd8"      "Tnfaip2"   "Marcksl1"  "Rep15"     "Cd79a"    
## [235] "Zfp503"    "Lig1"      "S100a10"   "Foxp3"     "Ccrl2"     "Tk1"      
## [241] "Serpinb2"  "Irf4"      "Klrg1"     "Csn3"      "Ms4a6b"    "Rgs6"     
## [247] "Cd40lg"    "Ccr1"      "Tsc22d1"   "Pmepa1"    "Bmp2"      "Gimap7"   
## [253] "Gbp2b"     "Dapl1"     "Tmem132c"  "Pbk"       "Npy"       "Tchh"     
## [259] "Pcdh18"    "Elavl2"    "Fam13a"    "Ezh2"      "Tmem158"   "Runx1t1"  
## [265] "Ar"        "Cxcr6"     "Cwc25"     "Ptma"      "Ccdc34"    "Klra9"    
## [271] "Eomes"     "Gzmk"      "Rgs10"     "Vgf"       "Serpina3f" "Adamts1"  
## [277] "Irf8"      "Rilpl2"    "Anp32b"    "Tagln"     "Rab27b"    "Tubb2b"   
## [283] "Lat"       "Anxa2"     "Ndrg1"     "Cd6"       "Gpx1"      "Nr4a1"    
## [289] "Got1"      "Lmnb1"     "Cidea"     "Epha5"     "Tnf"       "Klri2"    
## [295] "Gbp6"      "Cenpm"     "Olfr889"   "Hmgn3"     "Klf4"      "Sdc4"
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB) +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB)

well, there’s nearly no big difference in PC1/2, a little bit shift in PC3/4, comparing with initial run

DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)

decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)

PCs <- 1:27
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 10)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 2.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 22088
## Number of edges: 306841
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7802
## Number of communities: 31
## Elapsed time: 1 seconds

Run UMAP/tSNE

GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 998)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 15:27:31 UMAP embedding parameters a = 0.9922 b = 1.112
## 15:27:31 Read 22088 rows and found 27 numeric columns
## 15:27:31 Using Annoy for neighbor search, n_neighbors = 20
## 15:27:31 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 15:27:34 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\Rtmp2phpoV\file95fc4de13e2a
## 15:27:34 Searching Annoy index using 1 thread, search_k = 2000
## 15:27:39 Annoy recall = 100%
## 15:27:40 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 15:27:41 Initializing from normalized Laplacian + noise (using irlba)
## 15:27:41 Commencing optimization for 200 epochs, with 655162 positive edges
## 15:28:04 Optimization finished
DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno", label.size = 2.75) +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno", label.size = 2.75) +
  DimPlot(GEX.seur, reduction = "umap", label = T)

sort and test

GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
                                 levels = c(1,8,      # Tn
                                            27, 19,   # Treg
                                            16,       # Th17
                                            25,       # Th.CC
                                            17,       # Th1
                                            30,       # Th2
                                            29,       # ?
                                            7,        # iNKT
                                            21,       # MAIT
                                            20,5,     # gdT
                                            26,10,28, # ILC2
                                            24,       # ILC3.CC
                                            14,22,18,2,4,6,9,23,      # ILC3.Ncr1
                                            0,                        # ILC3.nn
                                            3,13,12,15,11             # ILC3.Ccr6
                                            ))
VlnPlot(subset(GEX.seur, subset = cnt %in% c("CTL")), 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "sort_clusters")

VlnPlot(subset(GEX.seur, subset = cnt %in% c("CKO")), 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "sort_clusters")

markers

# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "sort_clusters"

#GEX.markers.sort <- FindAllMarkers(GEX.seur, only.pos = TRUE, min.pct = 0.05,
#                                  test.use = "MAST",
#                                  #test.use = "wilcox",
#                                  logfc.threshold = 0.25)
GEX.markers.sort <- read.table("20230927_10x_SZJ.t1_adv.sort_markers.csv", header = TRUE, sep = ",")
#GEX.markers.sort %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
GEX.markers.sort$cluster <- factor(as.character(GEX.markers.sort$cluster),
                          levels = levels(GEX.seur$sort_clusters))



markers.sort_t60 <- (GEX.markers.sort %>% group_by(cluster) %>% 
                  filter(pct.1>0.01 & gene %in% grep("Rps|Rpl|mt-",GEX.markers.sort$gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                    filter(p_val_adj < 0.01) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.sort_t120 <- (GEX.markers.sort %>% group_by(cluster) %>% 
                  filter(pct.1>0.01 & gene %in% grep("Rps|Rpl|mt-",GEX.markers.sort$gene,invert = T,value = T)) %>%
                   top_n(n = 120, wt = avg_log2FC) %>%
                    filter(p_val_adj < 0.01) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
747/60
## [1] 12.45
747/64
## [1] 11.67188
747/65
## [1] 11.49231
pp.t60 <- list()

for(i in 1:12){
pp.t60[[i]] <- DotPlot(GEX.seur, features = rev(markers.sort_t60[(64*i-63):(64*i)]))  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}

pp.t60
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

nnn = 1278
nnn/60
## [1] 21.3
nnn/64
## [1] 19.96875
nnn/65
## [1] 19.66154
pp.t120 <- list()

for(i in 1:20){
pp.t120[[i]] <- DotPlot(GEX.seur, features = rev(markers.sort_t120[(64*i-63):(64*i)]))  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}
## Warning in FetchData.Seurat(object = object, vars = features, cells = cells):
## The following requested variables were not found: NA
#pp.t120
some reference papers for innate-like Tcell markers

Shilpi Chandra et al. Sci Immunol. 2023
https://pubmed.ncbi.nlm.nih.gov/37948512/

S Harsha Krovi NatCommun. 2020
https://pubmed.ncbi.nlm.nih.gov/33288744/

MAIT: Klrd1,Klra9,Nkg7,Ccl5,Cd7,Cd160,Xcl1 …
iNKT: Klrb1c,Ifng,Tbx21,Gzma …

still, here iNKT/MAIT clusters couldn’t be 100% confirmed,
because they usually are defined by specific TCR v-j genes

VlnPlot(GEX.seur, features = "Pdcd1", group.by = "sort_clusters") + NoLegend()

VlnPlot(GEX.seur, features = c("Fgl2","Klrd1","Klrk1","Klra9",
                               "Styk1","Nkg7","Ccl5","Cd7"), group.by = "sort_clusters",ncol = 1) + NoLegend()

VlnPlot(GEX.seur, features = c("Cd160","Klrb1c","Xcl1","Cxcr3","Gimap4"), group.by = "sort_clusters", ncol = 1) + NoLegend()

VlnPlot(GEX.seur, features = c("Klrb1c","Ifng","Tbx21","Gzma"), group.by = "sort_clusters",ncol = 1) + NoLegend()

VlnPlot(GEX.seur, features = "Eomes", group.by = "sort_clusters") + NoLegend()

DotPlot(GEX.seur, features = markers.new, group.by = "sort_clusters",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 9)) + 
  scale_y_discrete(limits=rev)

Anno

# remove very few low-quality mix-like C29    
GEX.seur <- subset(GEX.seur, subset = seurat_clusters != 29)
GEX.seur
## An object of class Seurat 
## 16621 features across 22038 samples within 1 assay 
## Active assay: RNA (16621 features, 1197 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
GEX.seur@meta.data[,grep("pANN|snn",colnames(GEX.seur@meta.data),value = T)] <- NULL
head(GEX.seur@meta.data)
##                    orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAACAGGA-1   ILC3_GEX       3483         1341   2.871088   29.08412
## AAACCCAAGAATTGCA-1   ILC3_GEX       2244         1068   4.099822   19.51872
## AAACCCAAGCAGCCCT-1   ILC3_GEX       2074          981   3.905497   23.09547
## AAACCCAAGCCTGAAG-1   ILC3_GEX       3232         1370   2.165842   22.92698
## AAACCCAAGCGCATCC-1   ILC3_GEX       2902         1288   2.446589   23.56995
## AAACCCACAGAAGTGC-1   ILC3_GEX       2306         1142   2.428448   21.50911
##                    FB.info      S.Score   G2M.Score Phase seurat_clusters cnt
## AAACCCAAGAACAGGA-1   CKO.2  0.009665415 -0.04022873     S              19 CKO
## AAACCCAAGAATTGCA-1   CTL.3 -0.050866345 -0.08275860    G1               6 CTL
## AAACCCAAGCAGCCCT-1   CKO.4 -0.069848995 -0.07870466    G1               0 CKO
## AAACCCAAGCCTGAAG-1   CKO.1 -0.040119167 -0.10621201    G1              20 CKO
## AAACCCAAGCGCATCC-1   CTL.2 -0.024552021 -0.03206281    G1              20 CTL
## AAACCCACAGAAGTGC-1   CTL.3 -0.068676919 -0.11434633    G1              13 CTL
##                    DoubletFinder0.05 DoubletFinder0.1 sort_clusters   preAnno
## AAACCCAAGAACAGGA-1           Singlet          Singlet            19      Treg
## AAACCCAAGAATTGCA-1           Singlet          Singlet             6 ILC3.Ncr1
## AAACCCAAGCAGCCCT-1           Singlet          Singlet             0   ILC3.nn
## AAACCCAAGCCTGAAG-1           Singlet          Singlet            20       gdT
## AAACCCAAGCGCATCC-1           Singlet          Singlet            20       gdT
## AAACCCACAGAAGTGC-1           Singlet          Singlet            13 ILC3.Ccr6
GEX.seur$Anno <- as.character(GEX.seur$seurat_clusters)

## Cd3d+Cd4+
GEX.seur$Anno[GEX.seur$Anno %in% c(1,8)] <- "Tn"     # Cd4+(weak)Sell+Ccr7+
GEX.seur$Anno[GEX.seur$Anno %in% c(27,19)] <- "Treg"   # Foxp3+
GEX.seur$Anno[GEX.seur$Anno %in% c(16)] <- "Th17"   # Il17a+
GEX.seur$Anno[GEX.seur$Anno %in% c(25)] <- "Th.CC"     # cycling, most should be Th17
GEX.seur$Anno[GEX.seur$Anno %in% c(17)] <- "Th1"       # Gzmk+
GEX.seur$Anno[GEX.seur$Anno %in% c(30)] <- "Th2"       # Cd3+Cd4+ and co-express many ILC2 markers

## Cd3d+Cd4-
GEX.seur$Anno[GEX.seur$Anno %in% c(7)] <- "iNKT"   # Cd160+Tbx21+Klrb1c+
GEX.seur$Anno[GEX.seur$Anno %in% c(21)] <- "MAIT"  # Cd160+Klrd1+Cd7+Itgae+
GEX.seur$Anno[GEX.seur$Anno %in% c(20,5)] <- "gdT"   # Tcrg-C1+Cd163l1+Trdv4+, C20 Il17a+

## Cd3d-
GEX.seur$Anno[GEX.seur$Anno %in% c(26,10,28)] <- "ILC2"     # Gata3+Calca+Il5+

# Rorc+
GEX.seur$Anno[GEX.seur$Anno %in% c(24)] <- "ILC3.CC"   # cycling, should be Ncr1+
GEX.seur$Anno[GEX.seur$Anno %in% c(14,22,18,2,4,6,9,23)] <- "ILC3.Ncr1"  # Ncr1+
GEX.seur$Anno[GEX.seur$Anno %in% c(0,3)] <- "ILC3.nn"    # Ncr1-Ccr6- double negative
GEX.seur$Anno[GEX.seur$Anno %in% c(13,12,15,11)] <- "ILC3.Ccr6"  # Ccr6+

GEX.seur$Anno <- factor(GEX.seur$Anno,
                        levels = c("Tn","Treg","Th17","Th.CC","Th1","Th2",
                                   "iNKT","MAIT","gdT",
                                   "ILC2",
                                   "ILC3.CC","ILC3.Ncr1","ILC3.nn","ILC3.Ccr6"))
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno", label.size = 2.75) +
  DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters")

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno", label.size = 2.75) +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

Idents(GEX.seur) <- "Anno"

GEX.seur$rep <- paste0("rep",
                       gsub("CTL|CKO","",as.character(GEX.seur$FB.info)))

#saveRDS(GEX.seur, "./Baf53cre-AhrCKO.SI_ILC3CD4T.Anno.rds")

final markers

## define feature colors
# Cell2020     
material.heat = function(n){
  colorRampPalette(
    c(
      "#283593", # indigo 800
      "#3F51B5", # indigo
      "#2196F3", # blue
      "#00BCD4", # cyan
      "#4CAF50", # green
      "#8BC34A", # light green
      "#CDDC39", # lime
      "#FFEB3B", # yellow
      "#FFC107", # amber
      "#FF9800", # orange
      "#FF5722"  # deep orange
    )
  )(n)
}
pp.Anno <- DotPlot(GEX.seur, features = markers.new, group.by = "Anno",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 9)) + 
  scale_color_gradientn(colours  = material.heat(100)) +
  scale_y_discrete(limits=rev)

pp.Anno

markers.final <- c("Cd3d","Cd4","Ccr7","Sell",
                   "Hopx","Foxp3","Il10",
                   "Il17a","Nebl","Top2a","Mki67",
                   "Gzma","Gzmk","Ccr5","Il12rb2","Ifng",
                   "Tbx21","Klrb1c","Xcl1","Nkg7","Plac8","Ccl5",
                   "Cd160","Klrd1","Cd7","Itgae",
                   "Tcrg-C1","Cd163l1","Blk",
                   "Gata3","Calca","Il5","Il13",
                   "Ncr1","Gzmb","Klrk1","Fcer1g",
                   "Il22","Ccr6","Rorc","Il17f")
pn.Anno.a <- DotPlot(GEX.seur, features = markers.final, group.by = "Anno",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 9)) + 
  scale_color_gradientn(colours  = material.heat(100)) +
  scale_y_discrete(limits=rev)

pn.Anno.a

pn.Anno.b <- DotPlot(GEX.seur, features = markers.final, group.by = "Anno",
        cols = c("midnightblue","darkorange1")) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 9)) + 
  scale_color_gradientn(colours  = material.heat(100)) #+
  #scale_y_discrete(limits=rev)

pn.Anno.b

DEGs

#
Idents(GEX.seur) <- "cnt"
all.clusters <- levels(GEX.seur$Anno)
all.clusters
##  [1] "Tn"        "Treg"      "Th17"      "Th.CC"     "Th1"       "Th2"      
##  [7] "iNKT"      "MAIT"      "gdT"       "ILC2"      "ILC3.CC"   "ILC3.Ncr1"
## [13] "ILC3.nn"   "ILC3.Ccr6"
#
list_new <- list()
list_new[["All"]] <- all.clusters

list_new[["CD4T"]] <- all.clusters[1:6]
list_new[["ILC3"]] <- all.clusters[11:14]


for(NN in all.clusters){
  list_new[[NN]] <- NN
}

names_new <- names(list_new)
list_new
## $All
##  [1] "Tn"        "Treg"      "Th17"      "Th.CC"     "Th1"       "Th2"      
##  [7] "iNKT"      "MAIT"      "gdT"       "ILC2"      "ILC3.CC"   "ILC3.Ncr1"
## [13] "ILC3.nn"   "ILC3.Ccr6"
## 
## $CD4T
## [1] "Tn"    "Treg"  "Th17"  "Th.CC" "Th1"   "Th2"  
## 
## $ILC3
## [1] "ILC3.CC"   "ILC3.Ncr1" "ILC3.nn"   "ILC3.Ccr6"
## 
## $Tn
## [1] "Tn"
## 
## $Treg
## [1] "Treg"
## 
## $Th17
## [1] "Th17"
## 
## $Th.CC
## [1] "Th.CC"
## 
## $Th1
## [1] "Th1"
## 
## $Th2
## [1] "Th2"
## 
## $iNKT
## [1] "iNKT"
## 
## $MAIT
## [1] "MAIT"
## 
## $gdT
## [1] "gdT"
## 
## $ILC2
## [1] "ILC2"
## 
## $ILC3.CC
## [1] "ILC3.CC"
## 
## $ILC3.Ncr1
## [1] "ILC3.Ncr1"
## 
## $ILC3.nn
## [1] "ILC3.nn"
## 
## $ILC3.Ccr6
## [1] "ILC3.Ccr6"
names_new
##  [1] "All"       "CD4T"      "ILC3"      "Tn"        "Treg"      "Th17"     
##  [7] "Th.CC"     "Th1"       "Th2"       "iNKT"      "MAIT"      "gdT"      
## [13] "ILC2"      "ILC3.CC"   "ILC3.Ncr1" "ILC3.nn"   "ILC3.Ccr6"

run

#test.DEGs_new
# save DEGs' table
df_test.DEGs_new <- data.frame()
for(nn in names_new){

  df_test.DEGs_new <- rbind(df_test.DEGs_new, data.frame(test.DEGs_new[[nn]],Anno=nn))
  
  
}

#write.csv(df_test.DEGs_new, "Baf53cre_AhrCKO.SI_ILC3CD4T.DEGs_CTLvsCKO_Anno.csv")

stat

cut0: padj 0.05 log2FC 0.1
cut1: padj 0.05 log2FC 0.2
cut2: padj 0.01 log2FC > log2(1.5) (|FC|>1.5)

cut0

# cut0
cut.padj = 0.05
cut.log2FC = 0.1  
  
cut.pct1 = 0.1
stat_test0.DEGs_new <- df_test.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(Anno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test0.DEGs_new
##         Anno cluster up.DEGs
## 1        All     CTL      50
## 2        All     CKO      19
## 3       CD4T     CTL       4
## 4       CD4T     CKO       8
## 5        gdT     CKO       4
## 6       ILC2     CKO       1
## 7       ILC3     CTL      39
## 8       ILC3     CKO      19
## 9  ILC3.Ccr6     CKO       3
## 10 ILC3.Ncr1     CTL      26
## 11 ILC3.Ncr1     CKO      19
## 12   ILC3.nn     CTL       3
## 13   ILC3.nn     CKO       3
## 14      iNKT     CKO       1
## 15       Th1     CTL       1
## 16        Tn     CTL       1
## 17        Tn     CKO       3
df_test.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(Anno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=Anno, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.cnt, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

cut1

# cut1
cut.padj = 0.05
cut.log2FC = 0.2  
  
cut.pct1 = 0.1
stat_test1.DEGs_new <- df_test.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(Anno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test1.DEGs_new
##         Anno cluster up.DEGs
## 1        All     CTL       4
## 2        All     CKO       3
## 3       CD4T     CTL       3
## 4       CD4T     CKO       6
## 5        gdT     CKO       2
## 6       ILC2     CKO       1
## 7       ILC3     CTL       6
## 8       ILC3     CKO       3
## 9  ILC3.Ccr6     CKO       2
## 10 ILC3.Ncr1     CTL       5
## 11 ILC3.Ncr1     CKO       5
## 12   ILC3.nn     CTL       3
## 13   ILC3.nn     CKO       1
## 14      iNKT     CKO       1
## 15       Th1     CTL       1
## 16        Tn     CTL       1
## 17        Tn     CKO       2
df_test.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(Anno,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=Anno, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.cnt, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

plot

pp_test.DEGs.new <- lapply(list_new, function(x){NA})
#
test.DEGs_new.combine <- test.DEGs_new
test.DEGs_new.combine <- lapply(test.DEGs_new.combine, function(x){
  x[x$cluster=="CTL","avg_log2FC"] <- -x[x$cluster=="CTL","avg_log2FC"]
  x
})

All

pp_test.DEGs.new$All <- test.DEGs_new.combine$All %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="All CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$All

CD4T

pp_test.DEGs.new$CD4T <- test.DEGs_new.combine$CD4T %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="CD4T CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$CD4T

ILC3

pp_test.DEGs.new$ILC3 <- test.DEGs_new.combine$ILC3 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="ILC3 CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC3

Tn

pp_test.DEGs.new$Tn <- test.DEGs_new.combine$Tn %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="Tn CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Tn

Treg

pp_test.DEGs.new$Treg <- test.DEGs_new.combine$Treg %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="Treg CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Treg

Th17

pp_test.DEGs.new$Th17 <- test.DEGs_new.combine$Th17 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="Th17 CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Th17

Th.CC

pp_test.DEGs.new$Th.CC <- test.DEGs_new.combine$Th.CC %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="Th.CC CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Th.CC

Th1

pp_test.DEGs.new$Th1 <- test.DEGs_new.combine$Th1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="Th1 CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Th1

Th2

pp_test.DEGs.new$Th2 <- test.DEGs_new.combine$Th2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="Th2 CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Th2

iNKT

pp_test.DEGs.new$iNKT <- test.DEGs_new.combine$iNKT %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="iNKT CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$iNKT

MAIT

pp_test.DEGs.new$MAIT <- test.DEGs_new.combine$MAIT %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="MAIT CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$MAIT

gdT

pp_test.DEGs.new$gdT <- test.DEGs_new.combine$gdT %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="gdT CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$gdT

ILC2

pp_test.DEGs.new$ILC2 <- test.DEGs_new.combine$ILC2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="ILC2 CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC2

ILC2

pp_test.DEGs.new$ILC2 <- test.DEGs_new.combine$ILC2 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="ILC2 CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC2

ILC3.CC

pp_test.DEGs.new$ILC3.CC <- test.DEGs_new.combine$ILC3.CC %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="ILC3.CC CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC3.CC

ILC3.Ncr1

pp_test.DEGs.new$ILC3.Ncr1 <- test.DEGs_new.combine$ILC3.Ncr1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="ILC3.Ncr1 CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC3.Ncr1

ILC3.nn

pp_test.DEGs.new$ILC3.nn <- test.DEGs_new.combine$ILC3.nn %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="ILC3.nn CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC3.nn

ILC3.Ccr6

pp_test.DEGs.new$ILC3.Ccr6 <- test.DEGs_new.combine$ILC3.Ccr6 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Gm",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
                               "CKO"=as.vector(color.cnt[2]),
                               "None"="grey")) +
  theme_classic() + labs(title="ILC3.Ccr6 CTL vs CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC3.Ccr6

composition

all

stat1_Anno <- GEX.seur@meta.data[,c("cnt","rep","Anno")]

stat1_Anno.s <- stat1_Anno %>%
  group_by(cnt,rep,Anno) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom1.Anno <- stat1_Anno.s %>%
  ggplot(aes(x = Anno, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.78, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.cnt, name="") +
  labs(title="Cell Composition of SI_ILC3CD4T.CTL_CKO, Anno", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=Anno, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom1.Anno

sig.test

glm.nb - anova.LRT

Sort.N <- c("CTL","CKO")

glm.nb_anova.1.Anno <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat1_Anno.s_N <- stat1_Anno.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat1_Anno.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat1_Anno.s_N$total <- total.N[paste0(stat1_Anno.s_N$cnt,
                                            "_",
                                            stat1_Anno.s_N$rep),"total"]
      
      glm.nb_anova.1.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat1_Anno.s_N$Anno)){
        glm.nb_anova.1.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat1_Anno.s_N %>% filter(Anno==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat1_Anno.s_N %>% filter(Anno==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat1_Anno.s_N %>% filter(Anno==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.1.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.1.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.1.Anno_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.1.Anno)))
rownames(glm.nb_anova.1.Anno_df) <- names(glm.nb_anova.1.Anno)
colnames(glm.nb_anova.1.Anno_df) <- gsub("X","C",colnames(glm.nb_anova.1.Anno_df))
glm.nb_anova.1.Anno_df
##                 Tn        Treg      Th17     Th.CC       Th1        Th2
## CTLvsCKO 0.3912709 0.008942821 0.8056086 0.8006229 0.7828195 0.03982599
##               iNKT      MAIT        gdT      ILC2    ILC3.CC ILC3.Ncr1
## CTLvsCKO 0.7431366 0.6684469 0.01670645 0.5941181 0.01892623 0.6151186
##            ILC3.nn ILC3.Ccr6
## CTLvsCKO 0.8142205 0.7319146
round(glm.nb_anova.1.Anno_df,4)
##              Tn   Treg   Th17  Th.CC    Th1    Th2   iNKT   MAIT    gdT   ILC2
## CTLvsCKO 0.3913 0.0089 0.8056 0.8006 0.7828 0.0398 0.7431 0.6684 0.0167 0.5941
##          ILC3.CC ILC3.Ncr1 ILC3.nn ILC3.Ccr6
## CTLvsCKO  0.0189    0.6151  0.8142    0.7319

CD4T

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno", label.size = 2.75, 
        cols = c(scales::hue_pal()(14)[c(1:6)],rep("lightgrey",8) )) +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

list_new$CD4T
## [1] "Tn"    "Treg"  "Th17"  "Th.CC" "Th1"   "Th2"
stat2_Anno <- GEX.seur@meta.data[,c("cnt","rep","Anno")] %>% 
                 filter(Anno %in% list_new$CD4T)

stat2_Anno.s <- stat2_Anno %>%
  group_by(cnt,rep,Anno) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom2.Anno <- stat2_Anno.s %>%
  ggplot(aes(x = Anno, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.78, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.cnt, name="") +
  labs(title="Cell Composition, CD4T.CTL_CKO, Anno", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=Anno, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom2.Anno

sig.test

glm.nb - anova.LRT

Sort.N <- c("CTL","CKO")

glm.nb_anova.2.Anno <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat2_Anno.s_N <- stat2_Anno.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat2_Anno.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat2_Anno.s_N$total <- total.N[paste0(stat2_Anno.s_N$cnt,
                                            "_",
                                            stat2_Anno.s_N$rep),"total"]
      
      glm.nb_anova.2.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat2_Anno.s_N$Anno)){
        glm.nb_anova.2.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat2_Anno.s_N %>% filter(Anno==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat2_Anno.s_N %>% filter(Anno==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat2_Anno.s_N %>% filter(Anno==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.2.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.2.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.2.Anno_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.2.Anno)))
rownames(glm.nb_anova.2.Anno_df) <- names(glm.nb_anova.2.Anno)
colnames(glm.nb_anova.2.Anno_df) <- gsub("X","C",colnames(glm.nb_anova.2.Anno_df))
glm.nb_anova.2.Anno_df
##                Tn         Treg      Th17     Th.CC       Th1         Th2 iNKT
## CTLvsCKO 0.736032 0.0006323642 0.9614593 0.3357168 0.1291926 0.002086836   NA
##          MAIT gdT ILC2 ILC3.CC ILC3.Ncr1 ILC3.nn ILC3.Ccr6
## CTLvsCKO   NA  NA   NA      NA        NA      NA        NA
round(glm.nb_anova.2.Anno_df,4)
##             Tn  Treg   Th17  Th.CC    Th1    Th2 iNKT MAIT gdT ILC2 ILC3.CC
## CTLvsCKO 0.736 6e-04 0.9615 0.3357 0.1292 0.0021   NA   NA  NA   NA      NA
##          ILC3.Ncr1 ILC3.nn ILC3.Ccr6
## CTLvsCKO        NA      NA        NA

ILCs

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno", label.size = 2.75, 
        cols = c(rep("lightgrey",9) ,scales::hue_pal()(14)[c(10:14)])) +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

list_new$ILC3
## [1] "ILC3.CC"   "ILC3.Ncr1" "ILC3.nn"   "ILC3.Ccr6"
stat3_Anno <- GEX.seur@meta.data[,c("cnt","rep","Anno")] %>% 
                 filter(Anno %in% c(list_new$ILC3,"ILC2"))

stat3_Anno.s <- stat3_Anno %>%
  group_by(cnt,rep,Anno) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom3.Anno <- stat3_Anno.s %>%
  ggplot(aes(x = Anno, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.78, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.cnt, name="") +
  labs(title="Cell Composition, ILCs.CTL_CKO, Anno", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=Anno, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom3.Anno

sig.test

glm.nb - anova.LRT

Sort.N <- c("CTL","CKO")

glm.nb_anova.3.Anno <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat3_Anno.s_N <- stat3_Anno.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat3_Anno.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat3_Anno.s_N$total <- total.N[paste0(stat3_Anno.s_N$cnt,
                                            "_",
                                            stat3_Anno.s_N$rep),"total"]
      
      glm.nb_anova.3.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat3_Anno.s_N$Anno)){
        glm.nb_anova.3.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat3_Anno.s_N %>% filter(Anno==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat3_Anno.s_N %>% filter(Anno==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat3_Anno.s_N %>% filter(Anno==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.3.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.3.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.3.Anno_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.3.Anno)))
rownames(glm.nb_anova.3.Anno_df) <- names(glm.nb_anova.3.Anno)
colnames(glm.nb_anova.3.Anno_df) <- gsub("X","C",colnames(glm.nb_anova.3.Anno_df))
glm.nb_anova.3.Anno_df
##          Tn Treg Th17 Th.CC Th1 Th2 iNKT MAIT gdT      ILC2    ILC3.CC
## CTLvsCKO NA   NA   NA    NA  NA  NA   NA   NA  NA 0.4784894 0.01651209
##          ILC3.Ncr1   ILC3.nn ILC3.Ccr6
## CTLvsCKO 0.3515031 0.6700376 0.7037525
round(glm.nb_anova.3.Anno_df,4)[,10:14,drop=F]
##            ILC2 ILC3.CC ILC3.Ncr1 ILC3.nn ILC3.Ccr6
## CTLvsCKO 0.4785  0.0165    0.3515    0.67    0.7038